

##function to get call rate###
args <- commandArgs(TRUE)
outpath<-args[2]
d<-args[1]
annfile<-args[3]
callrate <- function(x){ 
	#round(sum(x<0.05,na.rm=T)/sum(!is.na(x)),6)
	round(sum(x<0.05,na.rm=T)/length(x),6)
}

	samplefiles<-list.files(d)
#	samplefiles<-samplefiles[grep("BRCA",samplefiles)] ## skip any other files
#	samplefiles<-samplefiles[grep("-",samplefiles)]
	#samplefiles<-samplefiles
	#samplefiles<-samplefiles[-grep("DESCRIPTION.txt",samplefiles)]
	#samplefiles<-samplefiles[-grep("MANIFEST.txt",samplefiles)]
	#samplefiles<-samplefiles[-grep("README_DCC.txt",samplefiles)]
	samplenames<-c()
	for (i in 1:length(samplefiles)){
		samplefiles[i]<-paste(d,samplefiles[i],sep="/")
		
		if (i==1){
		
		#methy<-read.table(samplefiles[i],header = T, skip = 1, sep = "\t",stringsAsFactors = F)
		methy<-c()
		methy<-read.table(samplefiles[i],header = T, sep = "\t",stringsAsFactors = F)
	
		samplenames[i]<-colnames(methy)[2]
		colnames(methy)<-methy[1,]
		methy<-methy[2:nrow(methy),]
		if(length(methy$M_Number_Beads) ==0 && length(methy$M_STDERR) ==0 && length(methy$U_Number_Beads)==0 && length(methy$U_STDERR)==0)
		{
			methy$M_Number_Beads<-rep("NA",nrow(methy))
			methy$M_STDERR <-rep("NA",nrow(methy))
			methy$U_Number_Beads  <-rep("NA",nrow(methy))
			methy$U_STDERR   <-rep("NA",nrow(methy))
			methy<-methy[,c(1,2,9,10,3,11,12,4,5,6,7,8)]
		}
		#methy[,2]<-as.numeric(methy[,2])
		methy[,"Methylated_Signal_Intensity (M)"]<-as.numeric(methy[,"Methylated_Signal_Intensity (M)"])
		#methy[,5]<-as.numeric(methy[,5])
		methy[,"Un-Methylated_Signal_Intensity (U)"]<-as.numeric(methy[,"Un-Methylated_Signal_Intensity (U)"])
		#methy$beta<-round(methy[,2]/(methy[,2]+methy[,5]),4)
		methy$beta<-round(methy[,"Methylated_Signal_Intensity (M)"]/(methy[,"Methylated_Signal_Intensity (M)"]+methy[,"Un-Methylated_Signal_Intensity (U)"]),4)
		#methy<-methy[,c(1,13,2,5,12)]
		methy<-methy[,c("Composite Element REF","beta","Methylated_Signal_Intensity (M)","Un-Methylated_Signal_Intensity (U)","Detection_P_Value")]
		colnames(methy)<-paste(samplenames[i],c("CpG","Avg_Beta","Signal_B","Signal_A","Detection Pval"),sep=".")
		
		}else{
		tmp<-c()
		tmp<-read.table(samplefiles[i],header = T, sep = "\t",stringsAsFactors = F)
	
		samplenames[i]<-colnames(tmp)[2]
		colnames(tmp)<-tmp[1,]
		tmp<-tmp[2:nrow(tmp),]
		if(length(tmp$M_Number_Beads) ==0 && length(tmp$M_STDERR) ==0 && length(tmp$U_Number_Beads)==0 && length(tmp$U_STDERR)==0)
		{
			tmp$M_Number_Beads<-rep("NA",nrow(tmp))
			tmp$M_STDERR <-rep("NA",nrow(tmp))
			tmp$U_Number_Beads  <-rep("NA",nrow(tmp))
			tmp$U_STDERR   <-rep("NA",nrow(tmp))
			tmp<-tmp[,c(1,2,9,10,3,11,12,4,5,6,7,8)]
		}
		#tmp <- read.table(samplefiles[i],header = T, skip = 1, sep = "\t",stringsAsFactors = F)
		#tmp[,2]<-as.numeric(tmp[,2])
		#tmp[,5]<-as.numeric(tmp[,5])
		#tmp$beta<-round(tmp[,2]/(tmp[,2]+tmp[,5]),4)
		#tmp<-tmp[,c(1,13,2,5,12)]
		tmp[,"Methylated_Signal_Intensity (M)"]<-as.numeric(tmp[,"Methylated_Signal_Intensity (M)"])
		tmp[,"Un-Methylated_Signal_Intensity (U)"]<-as.numeric(tmp[,"Un-Methylated_Signal_Intensity (U)"])
		tmp$beta<-round(tmp[,"Methylated_Signal_Intensity (M)"]/(tmp[,"Methylated_Signal_Intensity (M)"]+tmp[,"Un-Methylated_Signal_Intensity (U)"]),4)
		tmp<-tmp[,c("Composite Element REF","beta","Methylated_Signal_Intensity (M)","Un-Methylated_Signal_Intensity (U)","Detection_P_Value")]
		colnames(tmp)<-paste(samplenames[i],c("CpG","Avg_Beta","Signal_B","Signal_A","Detection Pval"),sep=".")
		methy<-merge(methy,tmp,by.x=colnames(methy)[1],by.y=colnames(tmp)[1])
		}
	}

names(methy)[1]<-"CpGID" 
#average beta
avgbeta<-methy[,c(1,seq(2,dim(methy)[2],4))]
pvals<-methy[,c(1,seq(5,dim(methy)[2],4))]
#avg.pval<-apply(pvals[,-1],2,mean,na.rm=T)
#samcallrate<-apply(pvals[,-1],2,callrate)
pvals1<-as.data.frame(lapply(pvals[,2:ncol(pvals)],as.numeric))
avg.pval<-apply(pvals1,2,mean,na.rm=T)
samcallrate<-apply(pvals1,2,callrate)
sampledrop<-names(samcallrate[samcallrate<0.95])
#creating qc file
qc<-rep("Y",length(samplenames))
samplenames1<-samplenames
samplenames1<-gsub("\\.","-",samplenames)
qc<-cbind(samplenames1, qc)
if(length(sampledrop)>0)
{
	qc[which(samcallrate<0.95),2]<-"N"
}
colnames(qc)<-c("SAMPLE_NAME","QC_SAMPLE_FLAG")
write.table(qc,paste(outpath,"sample.qc",sep="/"),row.names=F,col.names=T,quote=F,sep="\t")

if(length(sampledrop)>0)
{
	sampledrop<-sub(".Detection.Pval","",sampledrop)
	sampledrop.3colms<-c(paste(sampledrop,".Avg_Beta",sep=""),paste(sampledrop,".Signal_B",sep=""),paste(sampledrop,".Signal_A",sep=""),paste(sampledrop,".Detection Pval",sep=""))	
	#only one sample would be dropped
	methy.keep<-methy[,-match(sampledrop.3colms,names(methy))]
}else{
	methy.keep<-methy
}
samcallrate.affilter<-apply(methy.keep[,c(1,seq(5,dim(methy.keep)[2],4))][,-1],2,callrate)
###read in the annotation file to get chr, CpG location#####
#methy27.anno<-read.csv("GPL8490_HumanMethylation27_270596_v.1.2.csv",skip = 7)
methy27.anno<-read.csv(annfile,skip = 7)
#get chr and location
methy27.anno.use<-methy27.anno[,match(c("IlmnID","Chr","MapInfo","Gene_Strand"),names(methy27.anno))]
names(methy27.anno.use)<-c("CpGID","Chr","Start","Strand")
methy27.anno.use$End<-methy27.anno.use$Start+1
methy27.anno.use$Strand <-paste(methy27.anno.use$Strand,"1",sep="")
#merge with methylation data###
methy.keep.anno<-merge(methy27.anno.use,methy.keep,by.x="CpGID",by.y="CpGID")

#write out one file for each sample###
#samnum<-(dim(methy.keep.anno)[2]-4)/4
samstart<-seq(6,dim(methy.keep.anno)[2],4)
kk_colnames<-c("CHR","MIN","MAX","STRAND","SEGID","METHYLATED","UNMETHYLATED","P_VALUE","BETA_VALUE")
i<-1
while(i<=length(samstart)){
	kk<-c()
	kk<-methy.keep.anno[,c(1:5,samstart[i]:(samstart[i]+3))]
	kk<-kk[,c(2,3,5,4,1,7,8,9,6)]
	#if(length(which(kk=="NA"))>1){kk<-gsub("NA","-",kk)}
	#if(length(which(kk=="+1"))>1){kk<-gsub("+1","1",kk)}
	#kk[kk == 'NA'] <- '-'
	kk[,6]<-round(kk[,6]);
	kk[,7]<-round(kk[,7]);
	kk[is.na(kk)]<-'-'
	kk[kk == '+1']<-'1'
	colnames(kk)<-kk_colnames
	kk$P_VALUE<-round(as.numeric(kk$P_VALUE),4)
	kk[is.na(kk)]<-'-'
	kk<-kk[with(kk, order(CHR)), ]
	write.table(kk,paste(outpath,paste(gsub("\\.","-",sub(".Avg_Beta","",colnames(methy.keep.anno)[samstart[i]])),"NORM",sep='.'),sep="/"),row.names=F,col.names=T,quote=F,sep="\t")
	#print(i)
	i<-i+1
}





